#define GPU inline __device__

// float functions
////////////////////////////////////////////////////////////////////////////////

GPU float lerp(float a, float b, float t)
{
    return a + t*(b-a);
}

GPU float clamp(float f, float a, float b)
{
    return fmaxf(a, fminf(f, b));
}

GPU float step(float edge, float x)
{
    return x < edge ? 0 : 1;
}

GPU float smoothstep(float edge0, float edge1, float x)
{
    float t = clamp((x - edge0) / (edge1 - edge0), 0, 1);
    return t * t * (3 - 2*t);
}

GPU float rect(float edge0, float edge1, float x)
{
    return edge0 <= x && x <= edge1 ? 1 : 0;
}



// float2 functions
////////////////////////////////////////////////////////////////////////////////

// additional constructors
GPU float2 make_float2(float s)
{
    return make_float2(s, s);
}
GPU float2 make_float2(int2 a)
{
    return make_float2(float(a.x), float(a.y));
}

// negate
GPU float2 operator-(float2 &a)
{
    return make_float2(-a.x, -a.y);
}

// addition
GPU float2 operator+(float2 a, float2 b)
{
    return make_float2(a.x + b.x, a.y + b.y);
}
GPU void operator+=(float2 &a, float2 b)
{
    a.x += b.x; a.y += b.y;
}

// subtract
GPU float2 operator-(float2 a, float2 b)
{
    return make_float2(a.x - b.x, a.y - b.y);
}
GPU void operator-=(float2 &a, float2 b)
{
    a.x -= b.x; a.y -= b.y;
}

// multiply
GPU float2 operator*(float2 a, float2 b)
{
    return make_float2(a.x * b.x, a.y * b.y);
}
GPU float2 operator*(float2 a, float s)
{
    return make_float2(a.x * s, a.y * s);
}
GPU float2 operator*(float s, float2 a)
{
    return make_float2(a.x * s, a.y * s);
}
GPU void operator*=(float2 &a, float s)
{
    a.x *= s; a.y *= s;
}

// divide
GPU float2 operator/(float2 a, float2 b)
{
    return make_float2(a.x / b.x, a.y / b.y);
}
GPU float2 operator/(float2 a, float s)
{
    float inv = 1.0f / s;
    return a * inv;
}
GPU float2 operator/(float s, float2 a)
{
    float inv = 1.0f / s;
    return a * inv;
}
GPU void operator/=(float2 &a, float s)
{
    float inv = 1.0f / s;
    a *= inv;
}

// lerp
GPU float2 lerp(float2 a, float2 b, float t)
{
    return a + t*(b-a);
}

// clamp
GPU float2 clamp(float2 v, float a, float b)
{
    return make_float2(clamp(v.x, a, b), clamp(v.y, a, b));
}

GPU float2 clamp(float2 v, float2 a, float2 b)
{
    return make_float2(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y));
}

// dot product
GPU float dot(float2 a, float2 b)
{ 
    return a.x * b.x + a.y * b.y;
}

// length
GPU float length(float2 v)
{
    return sqrtf(dot(v, v));
}

// normalize
GPU float2 normalize(float2 v)
{
    float invLen = rsqrtf(dot(v, v));
    return v * invLen;
}

// floor
GPU float2 floor(const float2 v)
{
    return make_float2(floor(v.x), floor(v.y));
}

// reflect
GPU float2 reflect(float2 i, float2 n)
{
	return i - 2.0f * n * dot(n,i);
}

// absolute value
GPU float2 fabs(float2 v)
{
	return make_float2(fabs(v.x), fabs(v.y));
}

// float3 functions
////////////////////////////////////////////////////////////////////////////////

// additional constructors
GPU float3 make_float3(float s)
{
    return make_float3(s, s, s);
}
GPU float3 make_float3(float2 a)
{
    return make_float3(a.x, a.y, 0.0f);
}
GPU float3 make_float3(float2 a, float s)
{
    return make_float3(a.x, a.y, s);
}
GPU float3 make_float3(float4 a)
{
    return make_float3(a.x, a.y, a.z);  // discards w
}
GPU float3 make_float3(int3 a)
{
    return make_float3(float(a.x), float(a.y), float(a.z));
}

// negate
GPU float3 operator-(float3 &a)
{
    return make_float3(-a.x, -a.y, -a.z);
}

// min
GPU float3 fminf(float3 a, float3 b)
{
	return make_float3(fminf(a.x,b.x), fminf(a.y,b.y), fminf(a.z,b.z));
}

// max
GPU float3 fmaxf(float3 a, float3 b)
{
	return make_float3(fmaxf(a.x,b.x), fmaxf(a.y,b.y), fmaxf(a.z,b.z));
}

// addition
GPU float3 operator+(float3 a, float3 b)
{
    return make_float3(a.x + b.x, a.y + b.y, a.z + b.z);
}
GPU float3 operator+(float3 a, float b)
{
    return make_float3(a.x + b, a.y + b, a.z + b);
}
GPU void operator+=(float3 &a, float3 b)
{
    a.x += b.x; a.y += b.y; a.z += b.z;
}

// subtract
GPU float3 operator-(float3 a, float3 b)
{
    return make_float3(a.x - b.x, a.y - b.y, a.z - b.z);
}
GPU float3 operator-(float3 a, float b)
{
    return make_float3(a.x - b, a.y - b, a.z - b);
}
GPU void operator-=(float3 &a, float3 b)
{
    a.x -= b.x; a.y -= b.y; a.z -= b.z;
}

// multiply
GPU float3 operator*(float3 a, float3 b)
{
    return make_float3(a.x * b.x, a.y * b.y, a.z * b.z);
}
GPU float3 operator*(float3 a, float s)
{
    return make_float3(a.x * s, a.y * s, a.z * s);
}
GPU float3 operator*(float s, float3 a)
{
    return make_float3(a.x * s, a.y * s, a.z * s);
}
GPU void operator*=(float3 &a, float s)
{
    a.x *= s; a.y *= s; a.z *= s;
}

// divide
GPU float3 operator/(float3 a, float3 b)
{
    return make_float3(a.x / b.x, a.y / b.y, a.z / b.z);
}
GPU float3 operator/(float3 a, float s)
{
    float inv = 1.0f / s;
    return a * inv;
}
GPU float3 operator/(float s, float3 a)
{
    float inv = 1.0f / s;
    return a * inv;
}
GPU void operator/=(float3 &a, float s)
{
    float inv = 1.0f / s;
    a *= inv;
}

// lerp
GPU float3 lerp(float3 a, float3 b, float t)
{
    return a + t*(b-a);
}

// clamp
GPU float3 clamp(float3 v, float a, float b)
{
    return make_float3(clamp(v.x, a, b), clamp(v.y, a, b), clamp(v.z, a, b));
}

GPU float3 clamp(float3 v, float3 a, float3 b)
{
    return make_float3(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y), clamp(v.z, a.z, b.z));
}

// dot product
GPU float dot(float3 a, float3 b)
{ 
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

// cross product
GPU float3 cross(float3 a, float3 b)
{ 
    return make_float3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); 
}

// length
GPU float length(float3 v)
{
    return sqrtf(dot(v, v));
}

// normalize
GPU float3 normalize(float3 v)
{
    float invLen = rsqrtf(dot(v, v));
    return v * invLen;
}

// floor
GPU float3 floor(const float3 v)
{
    return make_float3(floor(v.x), floor(v.y), floor(v.z));
}

// reflect
GPU float3 reflect(float3 i, float3 n)
{
	return i - 2.0f * n * dot(n,i);
}

// absolute value
GPU float3 fabs(float3 v)
{
	return make_float3(fabs(v.x), fabs(v.y), fabs(v.z));
}

// float4 functions
////////////////////////////////////////////////////////////////////////////////

// additional constructors
GPU float4 make_float4(float s)
{
    return make_float4(s, s, s, s);
}
GPU float4 make_float4(float a, float s)
{
    return make_float4(a, a, a, s);
}
GPU float4 make_float4(float3 a)
{
    return make_float4(a.x, a.y, a.z, 0.0f);
}
GPU float4 make_float4(float3 a, float w)
{
    return make_float4(a.x, a.y, a.z, w);
}
GPU float4 make_float4(int4 a)
{
    return make_float4(float(a.x), float(a.y), float(a.z), float(a.w));
}

// negate
GPU float4 operator-(float4 &a)
{
    return make_float4(-a.x, -a.y, -a.z, -a.w);
}

// min
GPU float4 fminf(float4 a, float4 b)
{
	return make_float4(fminf(a.x,b.x), fminf(a.y,b.y), fminf(a.z,b.z), fminf(a.w,b.w));
}

// max
GPU float4 fmaxf(float4 a, float4 b)
{
	return make_float4(fmaxf(a.x,b.x), fmaxf(a.y,b.y), fmaxf(a.z,b.z), fmaxf(a.w,b.w));
}

// addition
GPU float4 operator+(float4 a, float4 b)
{
    return make_float4(a.x + b.x, a.y + b.y, a.z + b.z,  a.w + b.w);
}
GPU void operator+=(float4 &a, float4 b)
{
    a.x += b.x; a.y += b.y; a.z += b.z; a.w += b.w;
}

// subtract
GPU float4 operator-(float4 a, float4 b)
{
    return make_float4(a.x - b.x, a.y - b.y, a.z - b.z,  a.w - b.w);
}
GPU void operator-=(float4 &a, float4 b)
{
    a.x -= b.x; a.y -= b.y; a.z -= b.z; a.w -= b.w;
}

// multiply
GPU float4 operator*(float4 a, float s)
{
    return make_float4(a.x * s, a.y * s, a.z * s, a.w * s);
}
GPU float4 operator*(float s, float4 a)
{
    return make_float4(a.x * s, a.y * s, a.z * s, a.w * s);
}
GPU void operator*=(float4 &a, float s)
{
    a.x *= s; a.y *= s; a.z *= s; a.w *= s;
}

// divide
GPU float4 operator/(float4 a, float4 b)
{
    return make_float4(a.x / b.x, a.y / b.y, a.z / b.z, a.w / b.w);
}
GPU float4 operator/(float4 a, float s)
{
    float inv = 1.0f / s;
    return a * inv;
}
GPU float4 operator/(float s, float4 a)
{
    float inv = 1.0f / s;
    return a * inv;
}
GPU void operator/=(float4 &a, float s)
{
    float inv = 1.0f / s;
    a *= inv;
}

// lerp
GPU float4 lerp(float4 a, float4 b, float t)
{
    return a + t*(b-a);
}

// clamp
GPU float4 clamp(float4 v, float a, float b)
{
    return make_float4(clamp(v.x, a, b), clamp(v.y, a, b), clamp(v.z, a, b), clamp(v.w, a, b));
}

GPU float4 clamp(float4 v, float4 a, float4 b)
{
    return make_float4(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y), clamp(v.z, a.z, b.z), clamp(v.w, a.w, b.w));
}

// dot product
GPU float dot(float4 a, float4 b)
{ 
    return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
}

// length
GPU float length(float4 r)
{
    return sqrtf(dot(r, r));
}

// normalize
GPU float4 normalize(float4 v)
{
    float invLen = rsqrtf(dot(v, v));
    return v * invLen;
}

// floor
GPU float4 floor(const float4 v)
{
    return make_float4(floor(v.x), floor(v.y), floor(v.z), floor(v.w));
}

// absolute value
GPU float4 fabs(float4 v)
{
	return make_float4(fabs(v.x), fabs(v.y), fabs(v.z), fabs(v.w));
}


// Frame
////////////////////////////////////////////////////////////////////////////////

class Frame
{
public:
    float3 x, y, z, origin;
    
    GPU void setDefault(float3 position)
    {
        origin = position;
        x = make_float3(1, 0, 0);
        y = make_float3(0, 1, 0);
        z = make_float3(0, 0, 1);
    }
    
    GPU void lookAt(float3 position, float3 target, float3 up)
    {
        origin = position;
        z = normalize(position - target);
        x = normalize(cross(up, z));
        y = normalize(cross(z, x));
    }
    
    GPU float3 getVectorToWorld(float3 v)
    {
        return v.x * x + v.y * y + v.z * z;
    }
    
    GPU float3 getPointToWorld(float3 p)
    {
        return p.x * x + p.y * y + p.z * z + origin;
    }
};

// transfer2
////////////////////////////////////////////////////////////////////////////////

GPU float transfer2(
    float x0, float f0,
    float x1, float f1,
    float x)
{
    if (x < x0) return 0;
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    return 0;
}

GPU float2 transfer2(
    float x0, float2 f0,
    float x1, float2 f1,
    float x)
{
    if (x < x0) return make_float2(0);
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    return make_float2(0);
}

GPU float3 transfer2(
    float x0, float3 f0,
    float x1, float3 f1,
    float x)
{
    if (x < x0) return make_float3(0);
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    return make_float3(0);
}

GPU float4 transfer2(
    float x0, float4 f0,
    float x1, float4 f1,
    float x)
{
    if (x < x0) return make_float4(0);
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    return make_float4(0);
}


// transfer3
////////////////////////////////////////////////////////////////////////////////

GPU float transfer3(
    float x0, float f0,
    float x1, float f1,
    float x2, float f2,
    float x)
{
    if (x < x0) return 0;
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    if (x < x2) return lerp(f1, f2, (x - x1) / (x2 - x0));
    return 0;
}

GPU float2 transfer3(
    float x0, float2 f0,
    float x1, float2 f1,
    float x2, float2 f2,
    float x)
{
    if (x < x0) return make_float2(0);
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    if (x < x2) return lerp(f1, f2, (x - x1) / (x2 - x0));
    return make_float2(0);
}

GPU float3 transfer3(
    float x0, float3 f0,
    float x1, float3 f1,
    float x2, float3 f2,
    float x)
{
    if (x < x0) return make_float3(0);
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    if (x < x2) return lerp(f1, f2, (x - x1) / (x2 - x0));
    return make_float3(0);
}

GPU float4 transfer3(
    float x0, float4 f0,
    float x1, float4 f1,
    float x2, float4 f2,
    float x)
{
    if (x < x0) return make_float4(0);
    if (x < x1) return lerp(f0, f1, (x - x0) / (x1 - x0));
    if (x < x2) return lerp(f1, f2, (x - x1) / (x2 - x0));
    return make_float4(0);
}


// helper textures for cubic interpolation and random numbers
////////////////////////////////////////////////////////////////////////////////

texture<float4, 2, cudaReadModeElementType> hgTexture;
texture<float4, 2, cudaReadModeElementType> dhgTexture;
texture<int, 2, cudaReadModeElementType> randomTexture;

GPU float3 hg(float a)
{
    // float a2 = a * a;
    // float a3 = a2 * a;
    // float w0 = (-a3 + 3*a2 - 3*a + 1) / 6;
    // float w1 = (3*a3 - 6*a2 + 4) / 6;
    // float w2 = (-3*a3 + 3*a2 + 3*a + 1) / 6;
    // float w3 = a3 / 6;
    // float g = w2 + w3;
    // float h0 = (1.0f + a) - w1 / (w0 + w1);
    // float h1 = (1.0f - a) + w3 / (w2 + w3);
    // return make_float3(h0, h1, g);
    return make_float3(tex2D(hgTexture, a, 0));
}

GPU float3 dhg(float a)
{
    return make_float3(tex2D(dhgTexture, a, 0));
}


// 3D data query functions
////////////////////////////////////////////////////////////////////////////////

#define INF            __int_as_float(0x7f800000)

GPU bool range_check(float3 p, float3 max)
{
    bool check1 = p.x >= 0 && p.y >= 0 && p.z >= 0;
    bool check2 = p.x <= max.x && p.y <= max.y && p.z <= max.z;
    return check1 && check2;
}

template<typename T> GPU float cubicQuery3DDirect(T data, float3 p)
{
    float3 alpha = 255 * (p - floor(p - 0.5f) - 0.5f);

    float3 hgx = hg(alpha.x);
    float3 hgy = hg(alpha.y);
    float3 hgz = hg(alpha.z);

    // 8 linear queries
    float q000 = tex3D(data, p.x - hgx.x, p.y - hgy.x, p.z - hgz.x);
    float q001 = tex3D(data, p.x - hgx.x, p.y - hgy.x, p.z + hgz.y);
    float q010 = tex3D(data, p.x - hgx.x, p.y + hgy.y, p.z - hgz.x);
    float q011 = tex3D(data, p.x - hgx.x, p.y + hgy.y, p.z + hgz.y);
    float q100 = tex3D(data, p.x + hgx.y, p.y - hgy.x, p.z - hgz.x);
    float q101 = tex3D(data, p.x + hgx.y, p.y - hgy.x, p.z + hgz.y);
    float q110 = tex3D(data, p.x + hgx.y, p.y + hgy.y, p.z - hgz.x);
    float q111 = tex3D(data, p.x + hgx.y, p.y + hgy.y, p.z + hgz.y);

    // lerp along z
    float q00 = lerp(q000, q001, hgz.z);
    float q01 = lerp(q010, q011, hgz.z);
    float q10 = lerp(q100, q101, hgz.z);
    float q11 = lerp(q110, q111, hgz.z);
    
    // lerp along y
    float q0 = lerp(q00, q01, hgy.z);
    float q1 = lerp(q10, q11, hgy.z);
    
    // lerp along x
    return lerp(q0, q1, hgx.z);
}

template<typename T> GPU float cubicQuery3DDirect(T data, float3 p, float3 size)
{
    if (!range_check(p, size)) return -INF;
    return cubicQuery3DDirect(data, p);
}

template<typename T> GPU float cubicQuery3D(T data, float3 p, float3 size)
{
    p = 0.5f * (p + 1.0f);
    p = p * size;
    return cubicQuery3DDirect(data, p, size);
}

template<typename T> GPU float linearQuery3DDirect(T data, float3 p)
{
    return tex3D(data, p.x, p.y, p.z);
}

template<typename T> GPU float linearQuery3DDirect(T data, float3 p, float3 size)
{
    if (!range_check(p, size)) return -INF;
    return tex3D(data, p.x, p.y, p.z);
}

template<typename T> GPU float linearQuery3D(T data, float3 p, float3 size)
{
    p = 0.5f * (p + 1.0f);
    p = p * size;
    return linearQuery3DDirect(data, p, size);
}


template<typename T> GPU float3 linearGradient3DDirect(T data, float3 p)
{
    float delta = 0.0001f;
    float dx = linearQuery3DDirect(data, make_float3(p.x + delta, p.y, p.z)) -
        linearQuery3DDirect(data, make_float3(p.x - delta, p.y, p.z));
    float dy = linearQuery3DDirect(data, make_float3(p.x, p.y + delta, p.z)) -
        linearQuery3DDirect(data, make_float3(p.x, p.y - delta, p.z));
    float dz = linearQuery3DDirect(data, make_float3(p.x, p.y, p.z + delta)) - 
        linearQuery3DDirect(data, make_float3(p.x, p.y, p.z - delta));
    return make_float3(dx, dy, dz) / (2 * delta);
}

template<typename T> GPU float3 linearGradient3DDirect(T data, float3 p, float3 size)
{
    if (!range_check(p, size)) return make_float3(0);
    return linearGradient3DDirect(data, p);
}

template<typename T> GPU float3 linearGradient3D(T data, float3 p, float3 size)
{
    if (!range_check(fabs(p), make_float3(1))) return make_float3(0);
    float delta = 0.01f;
    float dx = linearQuery3D(data, make_float3(p.x + delta, p.y, p.z), size) -
        linearQuery3D(data, make_float3(p.x - delta, p.y, p.z), size);
    float dy = linearQuery3D(data, make_float3(p.x, p.y + delta, p.z), size) -
        linearQuery3D(data, make_float3(p.x, p.y - delta, p.z), size);
    float dz = linearQuery3D(data, make_float3(p.x, p.y, p.z + delta), size) - 
        linearQuery3D(data, make_float3(p.x, p.y, p.z - delta), size);
    return make_float3(dx, dy, dz) / (2 * delta);
}


GPU float d_lerp(float a, float b, float t)
{
    return (b - a) * t;
}

template<typename T> GPU float3 cubicGradient3DDirect(T data, float3 p)
{
    float3 alpha = 255 * (p - floor(p - 0.5f) - 0.5f);

    float3 hgx = hg(alpha.x);
    float3 hgy = hg(alpha.y);
    float3 hgz = hg(alpha.z);

    float3 dhgx = dhg(alpha.x);
    float3 dhgy = dhg(alpha.y);
    float3 dhgz = dhg(alpha.z);

    // compute x-derivative
    float q000 = tex3D(data, p.x - dhgx.x, p.y - hgy.x, p.z - hgz.x);
    float q001 = tex3D(data, p.x - dhgx.x, p.y - hgy.x, p.z + hgz.y);
    float q010 = tex3D(data, p.x - dhgx.x, p.y + hgy.y, p.z - hgz.x);
    float q011 = tex3D(data, p.x - dhgx.x, p.y + hgy.y, p.z + hgz.y);
    float q100 = tex3D(data, p.x + dhgx.y, p.y - hgy.x, p.z - hgz.x);
    float q101 = tex3D(data, p.x + dhgx.y, p.y - hgy.x, p.z + hgz.y);
    float q110 = tex3D(data, p.x + dhgx.y, p.y + hgy.y, p.z - hgz.x);
    float q111 = tex3D(data, p.x + dhgx.y, p.y + hgy.y, p.z + hgz.y);

    float q00 = lerp(q000, q001, hgz.z);
    float q01 = lerp(q010, q011, hgz.z);
    float q10 = lerp(q100, q101, hgz.z);
    float q11 = lerp(q110, q111, hgz.z);
    
    float q0 = lerp(q00, q01, hgy.z);
    float q1 = lerp(q10, q11, hgy.z);
    
    float gradientX = d_lerp(q0, q1, dhgx.z);

    // compute y-derivative
    q000 = tex3D(data, p.x - hgx.x, p.y - dhgy.x, p.z - hgz.x);
    q001 = tex3D(data, p.x - hgx.x, p.y - dhgy.x, p.z + hgz.y);
    q010 = tex3D(data, p.x - hgx.x, p.y + dhgy.y, p.z - hgz.x);
    q011 = tex3D(data, p.x - hgx.x, p.y + dhgy.y, p.z + hgz.y);
    q100 = tex3D(data, p.x + hgx.y, p.y - dhgy.x, p.z - hgz.x);
    q101 = tex3D(data, p.x + hgx.y, p.y - dhgy.x, p.z + hgz.y);
    q110 = tex3D(data, p.x + hgx.y, p.y + dhgy.y, p.z - hgz.x);
    q111 = tex3D(data, p.x + hgx.y, p.y + dhgy.y, p.z + hgz.y);

    q00 = lerp(q000, q001, hgz.z);
    q01 = lerp(q010, q011, hgz.z);
    q10 = lerp(q100, q101, hgz.z);
    q11 = lerp(q110, q111, hgz.z);
    
    q0 = d_lerp(q00, q01, dhgy.z);
    q1 = d_lerp(q10, q11, dhgy.z);
    
    float gradientY = lerp(q0, q1, hgx.z);

    // compute z-derivative
    q000 = tex3D(data, p.x - hgx.x, p.y - hgy.x, p.z - dhgz.x);
    q001 = tex3D(data, p.x - hgx.x, p.y - hgy.x, p.z + dhgz.y);
    q010 = tex3D(data, p.x - hgx.x, p.y + hgy.y, p.z - dhgz.x);
    q011 = tex3D(data, p.x - hgx.x, p.y + hgy.y, p.z + dhgz.y);
    q100 = tex3D(data, p.x + hgx.y, p.y - hgy.x, p.z - dhgz.x);
    q101 = tex3D(data, p.x + hgx.y, p.y - hgy.x, p.z + dhgz.y);
    q110 = tex3D(data, p.x + hgx.y, p.y + hgy.y, p.z - dhgz.x);
    q111 = tex3D(data, p.x + hgx.y, p.y + hgy.y, p.z + dhgz.y);

    q00 = d_lerp(q000, q001, dhgz.z);
    q01 = d_lerp(q010, q011, dhgz.z);
    q10 = d_lerp(q100, q101, dhgz.z);
    q11 = d_lerp(q110, q111, dhgz.z);
    
    q0 = lerp(q00, q01, hgy.z);
    q1 = lerp(q10, q11, hgy.z);
    
    float gradientZ = lerp(q0, q1, hgx.z);

    return make_float3(gradientX, gradientY, gradientZ);
}

template<typename T> GPU float3 cubicGradient3DDirect(T data, float3 p, float3 size)
{
    if (!range_check(p, size)) return make_float3(0);
    return cubicGradient3DDirect(data, p);
}

template<typename T> GPU float3 cubicGradient3D(T data, float3 p, float3 size)
{
    p = 0.5f * (p + 1.0f);
    p = p * size;
    return cubicGradient3DDirect(data, p, size) * size;
}

// GPU float3 cubicGradient3D(T data, float3 p, float3 size)
// {
    // float delta = 0.01f;
    // float dx = cubicQuery3D(data, make_float3(p.x + delta, p.y, p.z), size) -
        // cubicQuery3D(data, make_float3(p.x - delta, p.y, p.z), size);
    // float dy = cubicQuery3D(data, make_float3(p.x, p.y + delta, p.z), size) -
        // cubicQuery3D(data, make_float3(p.x, p.y - delta, p.z), size);
    // float dz = cubicQuery3D(data, make_float3(p.x, p.y, p.z + delta), size) - 
        // cubicQuery3D(data, make_float3(p.x, p.y, p.z - delta), size);
    // return make_float3(dx, dy, dz) / (2 * delta);
// }


// 4D data query functions
////////////////////////////////////////////////////////////////////////////////

template<typename T> GPU float linearQuery4D(T data, float3 p, float4 size, float time)
{
    p = 0.5f * (p + 1.0f);
    p = p * make_float3(size);
    p.z += size.z * (int) time;
    return linearQuery3DDirect(data, p);
}

template<typename T> GPU float cubicQuery4D(T data, float3 p, float4 size, float time)
{
    p = 0.5f * (p + 1.0f);
    p = p * make_float3(size);
    p.z += size.z * (int) time;
    return cubicQuery3DDirect(data, p);
}

template<typename T> GPU float3 linearGradient4D(T data, float3 p, float4 size, float time)
{
    p = 0.5f * (p + 1.0f);
    p = p * make_float3(size);
    p.z += size.z * (int) time;
    return cubicGradient3DDirect(data, p) * make_float3(size);
}

template<typename T> GPU float3 cubicGradient4D(T data, float3 p, float4 size, float time)
{
    p = 0.5f * (p + 1.0f);
    p = p * make_float3(size);
    p.z += size.z * (int) time;
    return cubicGradient3DDirect(data, p) * make_float3(size);
}


// cut plane
////////////////////////////////////////////////////////////////////////////////

GPU float3 applyCutPlane(float3 rayStart, float3 rayEnd, float cut)
{
    float rayLength = length(rayEnd - rayStart);
    float3 omega = (rayEnd - rayStart) / rayLength;
    float t = cut - dot(rayStart, omega);
    t = clamp(t, 0, rayLength);
    return rayStart + t * omega;
}

GPU float4 makePlaneEquation(float3 p, float3 n)
{
    return make_float4(n, -dot(p, n));
}

GPU float evalPlaneEquation(float3 p, float4 eq)
{
    return p.x * eq.x + p.y * eq.y + p.z * eq.z + eq.w;
}


// Cut planes are specified by their plane equation (a float4).
// The normal of the plane points *into* the visible halfspace.

template<typename T> GPU float linearQuery3DCut(T data, float3 p, float3 size, float4 plane)
{
    if (evalPlaneEquation(p, plane) < 0) return -INF;
    return linearQuery3D(data, p, size);
}

template<typename T> GPU float3 linearGradient3DCut(T data, float3 p, float3 size,
    float4 plane, float step = 0.01f)
{
    float v = evalPlaneEquation(p, plane);
    if (v < 0) return make_float3(plane);
    float3 grad = linearGradient3D(data, p, size);
    if (v > step) return grad;
    return lerp(make_float3(plane), grad, v / step);
}

template<typename T> GPU float cubicQuery3DCut(T data, float3 p, float3 size, float4 plane)
{
    if (evalPlaneEquation(p, plane) < 0) return -INF;
    return cubicQuery3D(data, p, size);
}

template<typename T> GPU float3 cubicGradient3DCut(T data, float3 p, float3 size,
    float4 plane, float step = 0.01f)
{
    float v = evalPlaneEquation(p, plane);
    if (v < 0) return make_float3(plane);
    float3 grad = cubicGradient3D(data, p, size);
    if (v > step) return grad;
    return lerp(make_float3(plane), grad, v / step);
}

template<typename T> GPU float cubicQuery4DCut(T data, float3 p, float4 size, float4 plane,
    float time)
{
    if (evalPlaneEquation(p, plane) < 0) return -INF;
    return cubicQuery4D(data, p, size, time);
}

template<typename T> GPU float3 cubicGradient4DCut(T data, float3 p, float4 size,
    float4 plane, float time, float step = 0.01f)
{
    float v = evalPlaneEquation(p, plane);
    if (v < 0) return make_float3(plane);
    float3 grad = cubicGradient4D(data, p, size, time);
    if (v > step) return grad;
    return lerp(make_float3(plane), grad, v / step);
}

#define PARAM(param, value)
#define PARAM_RANGE(param, value, min, max, key)
#define DATA(name, file, type, dim, channels)


// definitions used in Python-generated code

#define linear_query_3d(d, x) linearQuery3D(d, x, d##_size)
#define linear_gradient_3d(d, x) linearGradient3D(d, x, d##_size)

#define linear_query_3d_direct(d, x) linearQuery3DDirect(d, x, d##_size)
#define linear_gradient_3d_direct(d, x) linearGradient3DDirect(d, x, d##_size)

#define linear_query_3d_cut(d, x, n, c) linearQuery3DCut(d, x, d##_size, make_float4(n, -c))
#define linear_gradient_3d_cut(d, x, n, c) linearGradient3DCut(d, x, d##_size, make_float4(n, -c))

#define linear_query_4d(d, x, t) linearQuery4D(d, x, d##_size, t)
#define linear_gradient_4d(d, x, t) linearGradient4D(d, x, d##_size, t)

#define linear_query_4d_cut(d, x, n, c, t) linearQuery4DCut(d, x, d##_size, make_float4(n, -c), t)
#define linear_gradient_4d_cut(d, x, n, c, t) linearGradient4DCut(d, x, d##_size, make_float4(n, -c), t)

#define cubic_query_3d(d, x) cubicQuery3D(d, x, d##_size)
#define cubic_gradient_3d(d, x) cubicGradient3D(d, x, d##_size)

#define cubic_query_3d_direct(d, x) cubicQuery3DDirect(d, x, d##_size)
#define cubic_gradient_3d_direct(d, x) cubicGradient3DDirect(d, x, d##_size)

#define cubic_query_3d_cut(d, x, n, c) cubicQuery3DCut(d, x, d##_size, make_float4(n, -c))
#define cubic_gradient_3d_cut(d, x, n, c) cubicGradient3DCut(d, x, d##_size, make_float4(n, -c))

#define cubic_query_4d(d, x, t) cubicQuery4D(d, x, d##_size, t)
#define cubic_gradient_4d(d, x, t) cubicGradient4D(d, x, d##_size, t)

#define cubic_query_4d_cut(d, x, n, c, t) cubicQuery4DCut(d, x, d##_size, make_float4(n, -c), t)
#define cubic_gradient_4d_cut(d, x, n, c, t) cubicGradient4DCut(d, x, d##_size, make_float4(n, -c), t)


#define linear_query_2d_rgba(d, p) tex2D(d, 0.5f * (p.x + 1.0f) * d##_size.x, 0.5f * (p.y + 1.0f) * d##_size.y)
#define linear_query_2d(d, p) tex2D(d, 0.5f * (p.x + 1.0f) * d##_size.x, 0.5f * (p.y + 1.0f) * d##_size.y)

#define linear_query_1d_rgba(d, x) tex2D(d, (0.5f * (x + 1.0f)) * d##_size.y, 0)
#define linear_query_1d(d, x) tex2D(d, (0.5f * (x + 1.0f)) * d##_size.y, 0)

#define pow powf
#define sin sinf
#define cos cosf
#define exp expf
#define log logf
#define T time


GPU float3 phong(float3 L, float3 N, float3 omega, float3 kd, float3 ks, float n, float amb)
{
    // diffuse
    float lobe = dot(N, L);
    if (lobe <= 0) return amb * kd;
    
    float3 color = (lobe + amb) * kd;
    
    // phong
    if (n > 0)
    {
        float3 R = reflect(-L, N);
        lobe = pow(fmaxf(dot(R, omega), 0), n);
        color += ks * lobe;
    }
    
    return color;
}

GPU float3 get_xyz(float4 v)
{
    return make_float3(v.x, v.y, v.z);
}






